#!/usr/bin/perl

##
## Copyright (C) 2008 Luciano Xumerle <luciano.xumerle@gmail.com>
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.
##

use strict;
use warnings;

#
# the genome must be in bfa format:
# use: maq fasta2bfa ref.fasta ref.bfa
# to convert your fasta sequence
# tricks: cat *fa | maq fasta2bfa - ref.bfa
#

my $maq = '/proj/genome/prog/maq-0.7.1/maq';
$maq = '/home/luciano/src/maq-0.7.1/maq'
  if ( -f '/home/luciano/src/maq-0.7.1/maq' );

my $gunzip = '/bin/gunzip';
$gunzip = '/usr/bin/gunzip' if ( -f '/usr/bin/gunzip' );

my $bunzip2 = '/bin/bunzip2';
$bunzip2 = '/usr/bin/bunzip2' if ( -f '/usr/bin/bunzip2' );

my $paired = 0;
my $genome = '';
my @solexa = ();
foreach (@ARGV)
{
    $genome = $_ if (m/\.bfa$/);

    $paired = 1 if (m/^-p$/);

    push @solexa, $_ if ( m/\.txt$/ || m/\.txt\.gz$/ || m/\.txt\.bz2$/ );
}

if ( $#solexa > -1 && $genome ne '' && $paired == 0 )
{
    foreach my $ff (@solexa)
    {
        &executeMAQ( $genome, $ff, '' );
    }
}
elsif ( $#solexa == 1 && $genome ne '' && $paired == 1 )
{
    &executeMAQ( $genome, $solexa[0], $solexa[1] );
}
else
{
    print q|runMaq.pl by Luciano Xumerle
SINTAX:
> runMaq.pl reference.bfa solexaFile1.txt [solexaFile2.txt, ...]
OR (paired-end)
> runMaq.pl reference.bfa -p solexaFile1.txt solexaFile2.txt
|;
}

#
# execute the maq pipeline
#
sub executeMAQ
{
    my $genome = shift;
    my $ori    = shift;
    my $ori2   = shift;

    my $solexa = &sol2bfq($ori);

    # Align the reads to the reference
    if ( $ori2 ne '' && -f $ori2 )
    {
        my $solexa2 = &sol2bfq($ori2);
        system( $maq, 'match', $solexa . '.map',
            $genome,
            $solexa . '.bfq',
            $solexa2 . '.bfq'
        );
    }
    else
    {
        system( $maq, 'match', $solexa . '.map', $genome, $solexa . '.bfq' );
    }

    # Statistics from the alignment
    `$maq mapcheck $genome $solexa.map > $solexa.mapcheck.txt`;

    # Build the mapping assembly
`$maq assemble  $solexa.consensus.cns $genome $solexa.map 2> $solexa.assemble.log`;

    # Extract consensus sequences and qualities
    `$maq cns2fq $solexa.consensus.cns >  $solexa.cns.fq`;

    # Extract list of SNPs
    `$maq cns2snp $solexa.consensus.cns > $solexa.cns.snp`;

    # Extract the Alignments
    `$maq mapview $solexa.map > $solexa.mapview.maq`;
}

sub sol2bfq
{
    my $ori    = shift;
    my $solexa = $ori;
    $solexa =~ s/.txt$//;
    $solexa =~ s/.txt.bz2$//;
    $solexa =~ s/.txt.gz$//;

    # convert solexa file to fastq format
    if ( $ori =~ m/\.txt$/ )
    {
        system( $maq, 'sol2sanger', $ori, $solexa . '.fastq' );
    }
    else
    {
        if ( $ori =~ m/\.txt\.gz$/ )
        {
            `$gunzip -c $ori | $maq sol2sanger - $solexa.fastq`;
        }
        else
        {
            `$bunzip2 -c $ori | $maq sol2sanger - $solexa.fastq`;
        }
    }

    # Convert the reads to the binary fastq format
    system( $maq, 'fastq2bfq', $solexa . '.fastq', $solexa . '.bfq' );
    return $solexa;
}

